#### Institutionalizing Mutual Toleration
#### Simone Wegmann
#### EJPR 2025
########################################################################################################

#set working directory here
setwd("")

library(foreign)    
library(MASS)			  
library(sandwich)		
library(zoo)			  
library(lmtest)			
library(Formula)		
library(betareg)		
library(mfx)			  
library(arm)
library(blme)
library(gtools)
library(lme4)
library(mlmRev)
library(memisc)
library(stargazer)
library(calibrate)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(haven)

library(tidyverse)
library(hrbrthemes)
library(viridis)
library(DT)
library(plotly)


#Data
Vdem <- read_dta("PPOPVDem.dta")
VdemH <- read_dta("PPOPVDem_dropH.dta")
Polity <- read_dta("PPOPPolity.dta")




################################################################################
## MAIN DOCUMENT
################################################################################

# Figure 1. Policy-Making Power of Opposition Players and Democratic Experience
#############################
par(oma=c(2,2,0.5,0.5),mar=c(2,3.5,2,1),mfrow=c(1,1))
plot(c(0:200), type="n", ylab="", xlab="", main="", xlim=c(0,200), ylim=range(0,1), lwd=2, yaxt="n", cex.lab=1.2, cex.axis=1.2, cex.main=1.2, las=1)
axis(2,at=c(0,0.2,0.4,0.6,0.8,1),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.2, las=1)
mtext(text="Democratic experience in years",side=1,outer=TRUE,pch=1,cex=1.2)
mtext(text="PPOP",side=2,outer=TRUE,cex=1.2)
points(Vdem$vdem4_age,Vdem$finalppop, main="", col=rgb(32,32,32,50,maxColorValue=255), pch=20, cex=1.1)




#Figure 2. Policy-Making Power of Opposition Players and Democratic Decline
#############################
plot(c(-0.5:1.5), type="n", ylab="", xlab="", main="", xlim=c(0,1), ylim=range(-0.5,1.5), lwd=2, yaxt="n", cex.lab=1.2, cex.axis=1.2, cex.main=1.2, mgp=c(2,.5,0))
axis(2,at=c(0,1),labels=c("Stability","Decline"), cex.axis=1.2)
abline(h=0, col="grey", lty = 1)
abline(h=1, col="grey", lty = 1)
mtext(text="Policy-Making Power of Opposition Players",side=1,outer=TRUE,pch=1,cex=1.2)

points(VdemH$finalppop,jitter(VdemH$demdec_4, 0.0), main="", col=rgb(32,32,32,50,maxColorValue=255), pch=19, cex=1.1)

Hungary <- Vdem[ which(Vdem$country=='Hungary'), ]
points(Hungary$finalppop,jitter(Hungary$demdec_4, 0.0), main="", col=rgb(32,32,32,50,maxColorValue=255), pch=17, cex=1.1)
text(Hungary$finalppop, Hungary$demdec_4, labels=Hungary$country, cex= 0.7, pos=3)




#Table 1. Democratic Decline
#############################
#M1
M1<- glm(Vdemdecline_4 ~ finalppopr
         + as.factor(vdem_lagvdem_dem):finalppopr
         + as.factor(vdem_lagvdem_dem)
         + timesinceso
         + vdem4_age + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
         data=VdemH, family=binomial("logit"))

coeftest(M1, vcov=vcovCL(M1, cluster=VdemH$gwno))


#M2
M2 <- glm(Vdemdecline_4 ~ finalppopr 
          + as.factor(vdem_lagvdem_dem):finalppopr 
          + as.factor(vdem_lagvdem_dem) 
          + timesinceso 
          + as.factor(vdem_agerec10) + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
          data=VdemH, family=binomial("logit"))
coeftest(M2, vcov=vcovCL(M2, cluster=VdemH$gwno))




#Figure 3. Predicted Probabilities of Democratic Decline
#############################
summary(M2)
coef(M2)

m.beta <- coef(M2)
v.beta <- vcovCL(M2, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,10,0,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

attach(mtcars)
par(oma=c(2,2,0.5,0.5),mar=c(2,1.5,2,1),mfrow=c(1,1))
par(xpd=F)

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.5, cex.axis=1.0, cex.main=1.5, cex.lab=1.5,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("weak","0.2","0.4","0.6","0.8","strong"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")

mtext(text="Policy-Making Power of Opposition Players",side=1,line=0,outer=TRUE,pch=1,cex=1.2)
mtext(text="Pr (democratic decline)",side=2,line=0,outer=TRUE,pch=1,cex=1.2)






################################################################################
## ONLINE APPENDIX
################################################################################


#Figure A.1
#############################
attach(mtcars)
par(oma=c(2,2,0.5,0.5),mar=c(2,1.5,2,1),mfrow=c(4,1))
par(xpd=F)

plot(c(-0.5:1.5), type="n", ylab="", xlab="", main="Vdem electoral democracy 4-point scale", xlim=c(0,1), ylim=range(-0.5,1.5), lwd=2, yaxt="n", cex.lab=1.2, cex.axis=1.2, cex.main=1.2, mgp=c(2,.5,0))
axis(2,at=c(0,1),labels=c("Stability","Decline"), cex.axis=1.1)
abline(h=0, col="grey", lty = 1)
abline(h=1, col="grey", lty = 1)
points(Vdem$finalppop,Vdem$demdec_4, main="", col=rgb(32,32,32,50,maxColorValue=255), pch=19, cex=1.1)

plot(c(-0.5:1.5), type="n", ylab="", xlab="", main="VDem substantial and significant change on continuous scale (5 year)", xlim=c(0,1), ylim=range(-0.5,1.5), lwd=2, yaxt="n", cex.lab=1.2, cex.axis=1.2, cex.main=1.2, mgp=c(2,.5,0))
axis(2,at=c(0,1),labels=c("Stability","Decline"), cex.axis=1.1)
abline(h=0, col="grey", lty = 1)
abline(h=1, col="grey", lty = 1)
points(Vdem$finalppop,Vdem$autoc52_demt1, main="", col=rgb(32,32,32,50,maxColorValue=255), pch=19, cex=1.1)

plot(c(-0.5:1.5), type="n", ylab="", xlab="", main="VDem substantial and significant change on continuous scale (10 year)", xlim=c(0,1), ylim=range(-0.5,1.5), lwd=2, yaxt="n", cex.lab=1.2, cex.axis=1.2, cex.main=1.2, mgp=c(2,.5,0))
axis(2,at=c(0,1),labels=c("Stability","Decline"), cex.axis=1.1)
abline(h=0, col="grey", lty = 1)
abline(h=1, col="grey", lty = 1)
points(Vdem$finalppop,Vdem$polydiff22_demt1, main="", col=rgb(32,32,32,50,maxColorValue=255), pch=19, cex=1.1)

plot(c(-0.5:1.5), type="n", ylab="", xlab="", main="Polity IV", xlim=c(0,1), ylim=range(-0.5,1.5), lwd=2, yaxt="n", cex.lab=1.2, cex.axis=1.2, cex.main=1.2, mgp=c(2,.5,0))
axis(2,at=c(0,1),labels=c("Stability","Decline"), cex.axis=1.1)
abline(h=0, col="grey", lty = 1)
abline(h=1, col="grey", lty = 1)
points(Polity$finalppop,Polity$demdecline, main="", col=rgb(32,32,32,50,maxColorValue=255), pch=19, cex=1.1)

mtext(text="Policy-Making Power of Opposition Players",side=1,outer=TRUE,pch=1,cex=1)




#Figure A.2
#############################
summary(M1)
coef(M1)

attach(mtcars)
par(oma=c(2,2,0.5,0.5),mar=c(2,1.5,2,1),mfrow=c(1,6))
par(xpd=F)

#1 year
m.beta <- coef(M1)
v.beta <- vcovCL(M1, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,11,1,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="1 year", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.lab=1.0,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")

#10 years
m.beta <- coef(M1)
v.beta <- vcovCL(M1, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,11,10,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="10 year", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.lab=1.0,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")



#20 years
m.beta <- coef(M1)
v.beta <- vcovCL(M1, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,11,20,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="20 year", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.lab=1.0,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")


#30 years
m.beta <- coef(M1)
v.beta <- vcovCL(M1, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,11,30,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="30 year", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.lab=1.0,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")


#40 years
m.beta <- coef(M1)
v.beta <- vcovCL(M1, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,11,40,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="40 year", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.lab=1.0,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")




#50 years
m.beta <- coef(M1)
v.beta <- vcovCL(M1, cluster=VdemH$gwno)
BETA <- mvrnorm(1000, m.beta, v.beta)
X <- cbind(1,seq(0,1,length.out=100),1,11,50,0,0,1,0,seq(0,1,length.out=100)*1)
y.hat <- X %*% t(BETA)
p.hat <- 1/(1+exp(-y.hat))
p.hat
p.ci <- apply(p.hat,1,quantile,probs=c(0.025,0.5,0.975))
p.ci

plot(c(1:100), type="n", ylab="", xlab="", xaxt="n", main="50 year", xlim=c(0,101), ylim=range(-0.01,1), lwd=2, cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.lab=1.0,mgp=c(2,.5,0))
axis(1,at=c(0,20,40,60,80,100),labels=c("0","0.2","0.4","0.6","0.8","1"), cex.axis=1.0)
abline(h=0, col="grey", lty = 1)

for (j in 1:1000){
  points(c(1:100),(1/(1+exp(-y.hat)))[,j], col=rgb(32,32,32,5,maxColorValue=255), type="l")
}
points(c(1:100),p.ci[2,], type="l", pch=20, lwd=1, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[1,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")
points(c(1:100),p.ci[3,], type="l", pch=20, lwd=1, lty=2, col="grey1", bty="n", xaxt="n")


mtext(text="Policy-Making Power of Opposition Players",side=1,line=0,outer=TRUE,pch=1,cex=1.0)
mtext(text="Pr (democratic decline)",side=2,line=0,outer=TRUE,pch=1,cex=1.0)



#Table A.1. Descriptive Statistics
#############################
variables <- c("Vdemdecline_4", "finalppopr", "presidential", "PR", "unit_rec", "GDPmedian", "timesinceso", "vdem_agerec10", "vdem4_age")

data_complete <- VdemH[complete.cases(VdemH[variables]), variables]

descriptives <- data.frame(
  variable = variables,
  N = sapply(data_complete, length),
  Mean = sapply(data_complete[variables], mean, na.rm = TRUE),
  SD = sapply(data_complete[variables], sd, na.rm = TRUE),
  Min = sapply(data_complete[variables], min, na.rm = TRUE),
  Max = sapply(data_complete[variables], max, na.rm = TRUE),
  Missing = sapply(data_complete[variables], function(x) sum(is.na(x)))
)

print(descriptives)


#Table A.2. Sample
#############################
data_complete_sample <- Vdem[complete.cases(Vdem[variables]), "country"]
sample <- sort(unique(data_complete_sample))
print(sample, n=50)



#Table A.3. Democratic Decline - Robustness Checks
#############################

####Model 1 
#M1.1: as M1 but with Hungary included 
M1.1<- glm(Vdemdecline_4 ~ finalppopr
           + as.factor(vdem_lagvdem_dem):finalppopr
           + as.factor(vdem_lagvdem_dem)
           + timesinceso
           + vdem4_age + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
           data=Vdem, family=binomial("logit"))
coeftest(M1.1, vcov=vcovCL(M1.1, cluster=Vdem$gwno))

#M1.2: as M1 but with Lührmann et al. measure, no Hungary
M1.2<- glm(autoc52 ~ finalppopr
           + as.factor(vdem_lagvdem_dem):finalppopr
           + as.factor(vdem_lagvdem_dem)
           + timesinceso
           + vdem4_age + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
           data=VdemH, family=binomial("logit"))
coeftest(M1.2, vcov=vcovCL(M1.2, cluster=VdemH$gwno))

#M1.3: as M1 but with Lührmann et al. measure, with Hungary
M1.3<- glm(autoc52 ~ finalppopr
           + as.factor(vdem_lagvdem_dem):finalppopr
           + as.factor(vdem_lagvdem_dem)
           + timesinceso
           + vdem4_age + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
           data=Vdem, family=binomial("logit"))
coeftest(M1.3, vcov=vcovCL(M1.3, cluster=Vdem$gwno))

####Model 2 
#M2.1: as M2 but with Hungary included 
M2.1 <- glm(Vdemdecline_4 ~ finalppopr 
          + as.factor(vdem_lagvdem_dem):finalppopr 
          + as.factor(vdem_lagvdem_dem) 
          + timesinceso 
          + as.factor(vdem_agerec10) + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
          data=Vdem, family=binomial("logit"))
coeftest(M2.1, vcov=vcovCL(M2.1, cluster=Vdem$gwno))

#M2.2: as M2 but with Lührmann et al. measure, no Hungary
M2.2 <- glm(autoc52 ~ finalppopr 
          + as.factor(vdem_lagvdem_dem):finalppopr 
          + as.factor(vdem_lagvdem_dem) 
          + timesinceso 
          + as.factor(vdem_agerec10) + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
          data=VdemH, family=binomial("logit"))
coeftest(M2.2, vcov=vcovCL(M2.2, cluster=VdemH$gwno))

#M2.3: as M2 but with Lührmann et al. measure, with Hungary
M2.3 <- glm(autoc52 ~ finalppopr 
            + as.factor(vdem_lagvdem_dem):finalppopr 
            + as.factor(vdem_lagvdem_dem) 
            + timesinceso 
            + as.factor(vdem_agerec10) + as.factor(presidential) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
            data=Vdem, family=binomial("logit"))
coeftest(M2.3, vcov=vcovCL(M2.3, cluster=Vdem$gwno))





#Table A.4. Democratic Decline. Interaction Presidentialism
#############################

#M1 
M1pres<- glm(Vdemdecline_4 ~ finalppopr
             + as.factor(vdem_lagvdem_dem):finalppopr
             + as.factor(vdem_lagvdem_dem)
             + as.factor(presidential):finalppopr
             + as.factor(presidential)
             + timesinceso
             + vdem4_age + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
             data=VdemH, family=binomial("logit"))
coeftest(M1pres, vcov=vcovCL(M1pres, cluster=VdemH$gwno))



#M2 
M2pres <- glm(Vdemdecline_4 ~ finalppopr 
              + as.factor(vdem_lagvdem_dem):finalppopr 
              + as.factor(vdem_lagvdem_dem) 
              + as.factor(presidential):finalppopr
              + as.factor(presidential)
              + timesinceso 
              + as.factor(vdem_agerec10) + as.factor(unit_rec) + as.factor(PR) + as.factor(GDPmedian),
              data=VdemH, family=binomial("logit"))
coeftest(M2pres, vcov=vcovCL(M2pres, cluster=VdemH$gwno))


